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Central 


Objectives 

Improve  pixel-level  detection:  Reduce  probability  of  false  alarm  for  given  Pd 

-  Address  optically-thick  plumes:  Improve  accuracy  of  estimated  path 
integrated  concentration  (column  density,  CL) 

-  Compatible  with  real-time  processing 

•  Limitations  of  current  practice 

-  Matched-filter-based  detection  presumes  optically-thin  plume 

-  Other  approaches  require  prior  measurements  of  background  -  not 
compatible  with  on-the-move  detection 


Payoff:  Improve  detection  immediately  following  large-scale 
release,  low-lying  plumes;  improve  mass  estimate 


Problem  Formulation 
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•  Ensemble  of  measured  spectra 

•  Measured  spectra  are  nonlinear  functions  of  atmospheric  temperature, 
constituent  profiles,  background  characteristics,  etc. 

•  Desire  inverse  solution  to  radiative  transfer  equation  (RTE) 

•  Inverse  solution  is  mathematically  ill-posed  -  no  unique  solution  for  Rn 
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Relation  to  Atmospheric  Profile  Retrieval 


Physical  Sciences  Inc. 


Ground  (Blackbody) 


Layer  2  j  Layer  1  | Background 


i 


(Atmosphere)  !  (Plume)  I  (Non-blackbody, 

Spatially-varying) 


•  Stratified  atmosphere  model 

•  Profile  retrieval 

-  Many  stratifications 

-  Simple  background 

-  Apply  constraints  to  layer-to-layer  variation 

•  Plume  detection 

-  Simple  atmosphere 

-  Complicated  background 

Apply  constraints  to  background  characterization 
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Simplified  Radiative  Transfer  Model 
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•  Simplifying  assumptions: 

Homogeneous  atmosphere  between  sensor  and  vapor  cloud 
-  Cloud  is  at  air  temperature 

•  Compare  performance  of  non-linear  (exact)  RT  model  with 
linearized  approximation 


Adaptive  InfraRed  Imaging  Spectroradiometer  (AIRIS) 
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•  Imaging  Fabry-Perot  spectrometer 

-  Mirror  spacing  ~  X 

-  Staring  IR  FPA 

-  Band  sequential  data  acquisition 
Co-registration  of  narrowband  images 

-  Tune  time- 2  ms 


•  Selective  sampling  of 
wavelengths 

Acquire  imagery  only  at 
wavelengths  which  facilitate 
target  ID 

Minimize  data  volume 

•  Wide  field -of-view,  wide 
spectral  coverage 


TEP  Detection: 

Shortcoming  of  Thin  Plume  Approximation 
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•  Triethyl  phosphate  (TEP)  release 

•  Post-processing: 

-Non-linear  estimator  in  IDL 

-False  alarm  mitigation:  4  of  8 
spatial  filter 

-Bad  pixels  substituted 

•  Detection  key: 

-TEP  only 
-Yellow:  OD  ~  0 
-Red:  OD  >  1 


B /2/200  6  1 6:59 :D0 , 1 50  C  ubfi060B02_5000 


Agenda 


Physical  Sciences  Inc. 


VG1 0-076-8 


•  Introduction 

•  Nonlinear  estimation 

-  Algorithm  formulation 

-  Test  data 

-  Results 


Conclusions 


Optimal  Estimation:  Bayesian  Approach 
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•  Bayesian  posterior  pdf  for  model  parameter  values: 
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•  Maximum  likelihood  parameter  values  maximize  posterior: 

i  6  -  arg  max{/?(0|5c )}  i 
i  =  arg  min{-  In  p{o\^c  )}i 

•  Multi-variate  normal  pdf  for  deviation  between  model  and  measurement: 

-  in  P(x\e) = \ [x  -  f(e)J d  '[x-  f(e)\ + cAe 

D  =  diag  |(T,2  ]} 

•  Prior  pdf  for  model  parameter  values 

-  In  p(6)  =  \\0-dJ[ R\0  -0a\+ce 


Optimal  Estimation:  Signal  Model 
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Signal  model:  f(@)~Te  °xo  Te\° La 

/  t  \  V 

calculated  calculated  estimated  Blackbody 
spectrum  plume  trans.  “baseline”  atTatm 


Plume  transmission: 


reference  spectrum 


Infrared  background: 

-  Linear  mixing  model  x0=/u  +  B/3 

-  Probabilistic  Principal  Components 

Robust  estimate  of  sample  covariance  (Huber-type  M -estimator) 

Model  parameters: 


0  =  \a,  Ta,0\ 


-  a:  Plume  OD 

-  Ta:  Plume/air  temperature 
(3:  Parameters  which  account  for  bkgd.  radiance  given  bkgd.  model 


Minimize  Cost  Function 
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•  Maximum  likelihood  parameter  values  minimize  cost  function 

•  Multivariate  normal  pdfs  result  in  "quadratic"  cost  function 


deviation  between  measured  and  deviation  of  parameters 

model  spectra  from  nominal  values 


T 

-  Quadratic  formulation:  C  =  r  r 

Prior  applied  to  background  coefficients  only: 

-  Residuals  vector:  r  =  [d  1/2  [x  —  /(#)];/?] 


[0-0jR.[0-9a\  =  PTP 


•  Determine  maximum  likelihood  parameter  values  by  nonlinear  estimation 

-  Approach  not  limited  to  quadratic  cost  function 

-  Quadratic  cost  function  amenable  to  computationally-efficient  solution 


Nonlinear  Optimization  Algorithms 
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•  Iterative  determination  of  parameters,  e.g.,  Newton's  Method: 

cost  function 


9a,=e,-H:‘(vc); 


Model  parameters 


iteration  no.  Hessian  matrix  gradient  vector 


•  Gauss-Newton  algorithm 

-  Appoximate  Hessian  matrix: 

-  Parameter  update  equation: 

Initial  guess  at  0  from  linear  model 


H  *  2JrJ 
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Levenberg-Marquardt  algorithm  also  applicable 
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Next  generation  algorithm(s) 


Test  Regions 
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Region  1 
(2.6  ±  0.5  K) 


AIRIS-WAD  datacube:  256  x  256  pixels 
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•  Plume-free  data  augmented  with 
synthetic  plumes: 

-  64  x  5  pixels 

Max  OD  from  0  to  3.0  (base  e) 

-  T(plume)  =  T(air)  =  25.0  deg  C 


•  Thermal  contrast 

- OK  along  horizon 

Monotonic  increase  with  elev.  angle 

•  Test  both  favorable  and 
unfavorable  regions 


Simulation:  Synthetic  R-134a  Plumes 
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•  Effective  plume  transmission: 

Reference  spectrum  from  PNNL 
library 

r  W  -  exp[-  CL  •  cr(/l)]i 

-  Specify  column  density 

Beer's  Law  +  instrument  resolution 
function 

•  Data  augmentation: 

-  Partition  measurement  into 
estimated  signal,  noise 

Modify  signal  w/plume  signature 

-  Add  back  estimated  noise 

i*P  =*o  +  [ i-7J°[4-*o]+£; 


i  = |  W  -gi^C-d' i 


Probability  of  Exceedance 
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Performance  Metric:  ROC  Curves 
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•  Binary  decision  hypotheses 

H0  ("plume  absent")  and  H1  ("plume  present") 

-  pdfs  for  detection  statistic:  p(F  \  //) 

•  ROC  curve  is  Pd(Fth)  vs  Pfa(Fth) 

-  Pd  from  plume-augmented  region 

-  Pfa  from  rest  of  scene 

•  ROC  "surface":  Pd(a;  Fth) 


uu 

PfaiF,h)=  \p(F\Ht))l^ 


F,u 


P„(F,h)=  \p(F\H,yin 


F,i 
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Performance  Comparison  with  Matched  Filter 
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•  Objective:  Compare  nonlinear  estimation  with  matched  filter  estimation 

-  Detection  statistics 

-  Column  density/optical  density 

•  Detection  with  nonlinear  estimator:  F  test 


r 


M)- 


•  Analogous  metric  for  clutter-matched  filter:  Adaptive  Cosine  Estimator 
(ACE) 

D  fe)=,  >'^-‘[^-4 _ 

MF  (.v'7  E  V)([x  -juf  X  [x  -//]) 


•  Matched-filter  optical  density  estimate: 


aMF  ~ 


s'T  ±~l{x  -  ju)  A T0 


\T  v1  —  1  t 

s  Zj  s 


AT 


eff 


Expect  near  identical  results  for  optically-thin  plumes 
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R-134a  Detection:  Optically-Thin  Plume,  OD=0.1 
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Region  1 
(2.6  ±  0.5  K) 


AIRIS-WAD  datacube:  256  x  256  pixels 


•  Plume  column  density  =  82  mg/m2  (20  ppmv-m) 

•  Detection  statistics  not  favorable  in  either  Region 

ACE  and  Gauss-Newton  ROC  curves  are  nearly  identical 

-  20  bands  in  test  datacube 

-  OD=0  reference  spectrum 
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R-134a  Detection:  Optically-Thin  Plume,  OD=0.3 
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Region  1 
(2.6  ±  0.5  K) 


AIRIS-WAD  datacube:  256  x  256  pixels 


Plume  column  density  =  246  mg/m2  (59  ppmv-m) 

Detection  statistics  not  favorable  in  Region  1,  marginal  in  Region  2 

-  Lower  thermal  contrast 

-  ~2  orders  of  magnitude  reduction  in  Pfa  from  Region  1  to  Region  2 

ACE  and  Gauss-Newton  ROC  curves  are  nearly  identical 


R-134a  Detection:  Optically-Thick  Plume,  OD=1.0 
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Region  1 
(2.6  ±  0.5  K) 


AIRIS-WAD  datacube:  256  x  256  pixels 
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False  Alarm  Rate 


•  Plume  column  density  =  822  mg/m2  (197  ppmv-m) 

•  Detection  statistics  favorable  in  Region  2,  marginal  in  Region  1 

-  >2  orders  of  magnitude  reduction  in  Pfa  from  Region  1  to  Region  2 

•  Gauss-Newton  produces  significantly  more  favorable  ROC  curves 
than  ACE 


Factor  of  ~2  improvement  in  Region  1  (Pfa  for  fixed  Pd) 
Multiple  orders  of  magnitude  improvement  in  Region  2 


R-134a  Detection:  Optically-Thick  Plume,  OD=2.0 

Physical  Sciences  Inc. 


Region  1 
(2.6  ±  0.5  K) 


AIRIS-WAD  datacube:  256  x  256  pixels 
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•  Plume  column  density  =  1643  mg/m2  (394  ppmv-m) 

•  Detection  statistics  favorable  in  both  Regions 

•  Gauss-Newton  produces  significantly  more  favorable  ROC  curves  than  ACE 

-  >1  order  of  magnitude  improvement  in  Region  1 

-  Multiple  orders  of  magnitude  improvement  in  Region  2 


Estimated  Peak  Optical  Density 


Column  Density  Estimation 
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Peak  Optical  Density  Peak  Optical  Density 

•  Increased  thermal  contrast  reduces  uncertainty,  no  effect  overall  accuracy 

•  Nonlinear  estimation 

-  Accurately  recovers  embedded  OD  (CL) 

-  Systematic  deviation  at  OD>1  is  instrument  resolution  effect 

•  Matched  Filter  systematically  underestimates  CL 

•  Nonlinear  estimator  always  as  good  or  better  than  MF 


Algorithm  Execution 
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•  Gauss-Newton  algorithm  is 
iterative 

•  Termination  criterion: 


0< 


r 

C, 


< s : 


max 


•  Initial  guess  is  Iteration  0 


•  Typical  results: 

1-2  iterations  for  no  plume  (plume 
OD=0) 

-  3  iterations  to  converge  for 
OD-2-3  TEP  plume 
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Background 
Plume,  OD=0.0 
Plume,  OD=1.0 
Plume,  OD=2.0 
Plume,  OD=3.0 


No.  Iterations 


•  Decreasing  5  to  0.0001  increase 
no.  iteration  but  no  statistically- 
significant  effect  on  CL 


Summary  and  Conclusions 
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•  Developed  nonlinear  estimator  for  plume  detection  and  characterization 
based  on  RTE 


-  Bayesian  formulation 

-  Statistical  model  for  IR  background 

-  Gauss-Newton  algorithm  to  estimate  maximum  a  posteriori  (MAP)  values 

•  Signal  model  developed  for  non-scattering  atmosphere,  single  layer  plume 

-  Easily  modified  to  address  more  complicated  atmospheres 

•  Nonlinear  estimation  significantly  outperforms  matched-filter-based  with 
optically-thick  plumes 

-  "Orders  of  magnitude"  improvement 

NL  estimator  and  matched  filter  produce  equivalent  results  for  optically-thin  plumes 

•  This  work  was  performed  under  Contracts  from  the  Defense  Threat 
Reduction  Agency  (HDTRA01-07-C-0067)  and  US  Army  ECBC  Aberdeen 
Proving  Ground,  MD  (W911SR-06-C-0022).  Any  opinions,  findings  and 
conclusions  or  recommendations  expressed  in  this  material  are  those  of 
the  author  and  do  not  necessarily  reflect  the  views  of  HDRA 

or  the  Army. 


Additional  Material 


Data  Processing  Chain 
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J-6549a 


•  Focus  is  pixel-level  target  detection 

•  New  background  characterization  approach  facilitates 
improved  pixel-level  detection 

•  "A  chain  is  only  as  strong  as  its  weakest  link." 

Provide  higher  quality  input  to  False  Alarm  Mitigation  block 
-  False  Alarm  Mitigation  is  separate  issue 


Technical  Approach 
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•  Adapt  methodology  used  for  atmospheric  profile  retrieval  from 
space-based  sensor  data  (e.g.  AIRS,  IASI,  MODIS,  TES) 

-  Parameterize  Radiative  Transfer  Equation  (RTE) 

Apply  Estimation  Theory  to  determine  max.  likelihood  parameter  values 

-  Exploit  large  data  set:  utilize  ensemble  statistics 

•  Rationale: 

-  Physics-based  model  for  observations 

-  Statistically-justified  constraints 

Strong  theoretical  foundation  (see,  e.g.,  C.D. Rodgers,  Inverse  Methods  for 
Atmospheric  Sounding) 

•  Benefits 

-  Adaptable  framework 

-  Immediate  application  to  non-scattering  atmosphere 

-  Can  modify  RTE  to  address  more  complicated  atmospheres 


Linear  Models 
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Xp  =  X o  +  a' S' 


x0=jU  +  B/3  x0  =  vb  ~  N(ju,  S) 

"Structured  background"  "Unstructured  background" 

•  "Structured  Background" 

-  Values  of  (3  are  unconstrained 

-  Generalized  Likelihood  Ratio  Test: 


Typical  implementation:  B  =  eigenvectors  of  sample  covariance  matrix 


Pg  =I-b(btb)'bj 


•  "Unstructured  Background" 

-  vb  is  a  random  vector 
Adaptive  Cosine  Estimator: 


•  Survey  article:  Manolakis,  Marden,  &  Shaw,  "Hyperspectral  Image 
Processing  for  ATR  Applications,"  Lincoln  Lab  J.,  v.14  (2003) 


Pros  and  Cons  of  Linear  Approximation 
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•  Pro:  Matrix  multiplication  results  in  fast  computation 

-  All  spectra  in  ensemble  may  be  processed  in  parallel 

Major  computational  expense  is  diagonalization  of  sample  covariance  matrix 

-  AIRIS-WAD:  <150  ms  to  process  65536  twenty  element  spectra  for  four 
target  signatures  (using  2005  vintage  technology) 


•  Pro:  Detection  statistics  well-understood  for  Gaussian  noise 


•  Con:  Underlying  physical  assumptions  not  valid  for  detection 
scenarios  of  interest 

-  Mathematical  model  not  matched  to  physics 


Linear  approximation  to  Beer's  Law  can  introduce  significant  error 


Why  Gauss-Newton  Yields  Better  Results 
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Wavelength  [jim] 

Spectrum  augmented  with  OD=3.0  plume 
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Wavelength  [|im] 

Ratios  of  rms  residuals  in  plume  region,  OD=3.0  to  OD=0 


•  Model  is  matched  to  the  data 

•  Fit  residuals  are  systematically  larger  with  linear  model 

-  Result  of  least-squares  minimization 

Location  of  largest  residuals  highly  correlated  with  strongest  R-134a  absorption  features 


Adaptive  Infrared  Imaging  Spectroradiometer  - 
Wide  Area  Detector  (AIRIS-WAD) 
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•  Optical: 

-  256  x  256  pixels 

-  30  deg  x  30  deg  FOV 

-  spectral  coverage:  7.9  to 

1 1 .2  pm  at  ~0.1  pm  resolution 
(~1%oU) 

•  Datacubes: 

-  20  wavelengths 

user  selectable  Vs,  specified 
prior  to  mission 

•  Real-time  datacube 
processing:  up  to  3  Hz 

•  Detection  algorithm  history: 

-  GLRT:  Winter  2005-Spring 
2006 

-  ACE:  since  Spring  2006 


Detector  Cryocooler 


Tunable 

Etalon 


1-1421 


Calibration  Blackbody 


Hyperspectral  Background  Model 
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Probabilistic  Principal  Components- 
based 

-  M.E. Tipping  &  C.M. Bishop,  J.R. Statist. 
Soc.  B  (1999) 

Linear  mixing  model 

x  =  jU  +  B/3 

i 

i 

Eigenvalue-based  covariance  i 
regularization 

E  «  S  =  BBt  +  sD 
2  =  DV2(uAUt)dV2 
B  =  Dll2U  (A  -si  )1/2 
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Eigenvalue  Rank 


•  E  =  robust  estimate  of  sample 

covariance:  Huber-type  M-estimator 


Gauss-Newton  Algorithm 


Physical  Sciences  Inc. 


•  Follows  from  Newton's  method  -  simplifying  approximations 

•  Good  for  solving  weakly  nonlinear  equations 

•  Hessian  matrix: 


m 


21 

q= 1 


j 


m 


21 

<7=1 


^L^L  +  r 
80,  80k  q  80,80,. 

80,  80k 


=  2JtJ 

i  dr  \ 

i  J  - —  r 

■  dO 

Jacobian 


Gradient: 


VeC  =  2  JTr 
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Parameter  update  equation: 
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Initial  guess  at  0  from  linear  model 


